/V&29- C£' z ^ 


NASA Contractor Report 172241 

ICASE 


NASA-CR-1 72241 
19840002769 


SPECTRAL METHODS FOR TIME DEPENDENT PARTIAL 
DIFFERENTIAL EQUATIONS 


David Gottlieb 
and 

Eli Turkel 





To 



Contract Nos. NAS1-16394 and NAS1-17130 
October 1983 


INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING 
NASA Langley Research Center, Hampton, Virginia 23665 


Operated by the Universities Space Research Association 


NASA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23665 



i;ny i ,a iqqq 


f r r*\/ 

ur- 

,.. c 




SPECTRAL METHODS FOR TIME DEPENDENT PARTIAL DIFFERENTIAL EQUATIONS 


David Gottlieb* 

Institute for Computer Applications in Science and Engineering 

and 

Tel-Aviv University 
Eli Turkel** 

Institute for Computer Applications in Science and Engineering 

and 

Tel-Aviv University 


ABSTRACT 

The theory of spectral methods for time dependent partial differential 
equations is reviewed. When the domain is periodic Fourier methods are 
presented while for nonperiodic problems both Chebyshev and Legendre methods 
are discussed. The theory is presented for both hyperbolic and parabolic 
systems using both Galerkin and collocation procedures. While most of the 

review considers problems with constant coefficients the extension to 
nonlinear problems is also discussed. Some results for problems with shocks 
are presented. 
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INTRODUCTION 


We begin by describing how to construct spectral approximations to time- 
dependent mixed initial-boundary value problems. We shall study differential 
equations of the form 


3u 

at 


Lu + f 


V 0) = u o 


( 1 . 1 ) 


where for each t, u(t) belongs to a Hilbert space H such that u 
satisfies homogeneous boundary conditions. For simplicity we assume that L 
is an unbounded time independent linear operator. 

Numerical methods can be characterized by specifying a finite dimensional 
subspace B N C H and a projection operator : H -* B^. We require that 
the sequence {P^} satisfies 


lim !IP X7 u - u II 

XT N 

N-foo 


0 . 


We shall concentrate on semi-discrete approximations to (1.1), i.e., time 
is still a continuous variable. Such a semi-discrete approximation can be 
written as 


3t 


P N LP N U N + P N 


f 


( 1 - 2 ) 


U N (0) = P N U 0 


The numerical approximation (1.2) converges to the solution 


lim II u 
N-h» 


N 



u II 


0. 


where e B r 
of (1.1) if 


(1.3) 
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Combining (1.1) and (1.2) and assuming is independent of t, the error 

satisfies the equation 

lit “ P N U ) = P N LP N^°N " P N U ) + P N^ LP N “ L ) U * < 1 ' A ) 

Now, P N LP n is an operator from B N to B N and so can be viewed as a 
matrix. In particular, exp(P N LP N ) is well defined. Hence the solution to 
(1.4) can be written as 

t 

U N ” P N U = / exp ^ P N LP N (t ' T ^V LP N " L ) u ( T > dT - (1*5) 

We call a scheme consistent if 

lim IP (L - LP )ull = 0 0 < t < T, (1.6) 

N>oo 

while the scheme is stable if 

exp(P N LP n t) < K(t), (1.7) 

where K is independent of N, the dimension of B N , i.e., exp(P N LPjj t) is 
uniformly bounded for all 0 < t < T. It then follows from (1.5) that if a 
scheme is consistent and stable then the scheme converges. 

For spectral methods we choose B N as the finite space of polynomials 
(or trigonometric polynomials) of degree at most N. The rationale behind 
this choice is that one can approximate arbitrary functions f by such 
polynomials and the rate of convergence is only governed by the smoothness of 
the function f. Hence we hope to obtain highly accurate approximations to 
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the solutions of (1.1). Different choices of the projection operator Pfj 

lead to different subclasses of spectral methods. 

In these lectures we shall only consider one-dimensional differential 
equations . 


2. FOURIER METHODS 

We first consider periodic problems with period 2m. For this case it 
is natural to let Bjj be {e^ X } , -N < j < N. 

(a) Galerkin Method 

Let v(x) £ H then 


v(x) ='l a R e 

n =-°o ! 


inx 


( 2 . 1 ) 


The Galerkin method is characterized by the projection operator P^ where 


N 


V v = 7 a 

N XT 11 

n=-N 


inx 


( 2 . 2 ) 


We now rewrite the approximation (1.2) in the form 

Visit -%-£)-» 

or using the definition of P N 


8u 


f N T * inx^ 

- Lu N " f » e 3 - 


-N < n < N 


(2.3) 


where 
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N 


U N ‘ £ a „ (t)e 

n=-N 


inx 


U N (0) = P N U 0 ‘ 


This can be interpreted as a system of 2N+1 ordinary differential equations 
for the coefficients a^(t). Equivalently one can expand the solution to 
(1.1) in a finite Fourier series and then truncate Lu N . Hence, the Galerkin 
method is equivalent to solving the system (1.1) in Fourier space rather than 
physical space. 

An alternative basis is to expand u(x) in terms of cos nx, 

0 < n < N, and sin nx, 1 < n < N-l. This is equivalent to demanding that 

a = a in (2.1). 

-n n 

(b) Pseudospectral Method 

The pseudospectral or collocation method is defined by letting P^ 
be an interpolation operator. If f(x) is a periodic function then P N f is 
the trigonometric interpolation of f at the collocation points Xj, i.e., 

P N f(x p = f(x j ) and P N f E V 

The following sets of points are the most commonly used collocation points 

j = 0, • • • ,2N-1 (2.4a) 

y j 2N+1 


.1 = 0 , 


,2N. 


(2.4b) 
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The Xj are useful when operating with a FFT based on an even number of 
points while the yj are useful for an odd number of points. We shall only 
describe the collocation method based on the Xj. In this case the operator 
P N is given by 

2N-1 

P M f(x) = l f(x )g (x). (2.5) 

j=0 3 3 

The g.j(x) are trigonometric polynomials of degree at most N and 
These polynomials are given explicitly by 

X — X 

g,<*> “ ^ sin[N(x-x ,)] cot — 2 — ~ • (2.6) 


' V 


The fact that gj(x) is a trigonometric polynomial of degree N follows 
from the equivalent representation 


. N U(x-x.) 

ejW-S I )T e 3 

3 £=-N C £ 


(2.7) 


where c. = 1 (|£|*N), Cjj = c_jj = 2. Thus we can represent Pjjf(x) either as 


2N-1 


~ 4.11 4 

P„ f ( x) « -RT 7 l f(x )sin[N(x-x )]cot 
j=0 3 3 


X - X , 


( 2 . 8 ) 


using (2.5) or as 


2N-1 7 t, . 

P f(x) - l f(x ) w l — 

W j=0 3 £=-N C £ 


N , i£(x-x_.) 


N n <■ 2N-1 -Ux. 

V 1 i&x 1 V \ j 

£ ~ e 2N” ^ 

£=-N C £ j=0 3 


(2.9) 
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using (2.7). Defining 


. 2N-1 -iJtx . 

l - l *(*,)« J 


a £ 2 Nc ^ *' x -i 

* *. j=0 J 


( 2 . 10 ) 


(2.9) becomes 


P f = £ a^e 

n „ _ „ A 


( 2 . 11 ) 


When applying the pseudospectral Fourier method, either the explicit 

interpolatory formula (2.8) or the complex-Fourier representation (2.10) - 

(2.11) may be used. The operator L is a differential operator and so it is 

d k f(x ) 

useful to obtain in terras of f(x.). One way is simply to 

dx* 3 

differentiate (2.8) and to evaluate the resulting expression at the points 


d n u N (x ) 2N-1 

— zs— ' .L w 


— - (D u)., 

, n n j’ 

dx J 


( 2 . 12 ) 


where D n is an 2N x 2N matrix with elements 


(D -v ■ 


and u is the column vector 


u(x 0 ) 




Explicitly, 


(D i>jk - 


... x. - x. 

3+k 3 k 


V 2 (“l)* 5 cot 


.1 * k 


i = k 


(2.13) 
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(D 2>jk 


1/2 (-D 1+11+1 


2N 2 + 1 


sin 


2 *1 *k 


j * k 


j = k 


More generally 


D n = ( V 


(2.14) 


(2.15) 


which easily follows from the properties of gj(x). is a real, 

antisymmetric matrix. In general D 2 ^ is a real, symmetric matrix while 

°2k+l a rea l> antisymmetric matrix. 

Computationally, the evaluation of derivatives using (2.13) - (2.15) 
involves the multiplication of an 2N-component vector u by an 2N x 2N 
matrix, D n , which typically requires 0(N ) arithmetic operations. However, 
since the matrix product is actually a convolutional sum, it is possible to 
use the FFT to evaluate (2.13) - (2.15) in only order N log N operations 
when N is a highly composite integer (like 2 P or 3^). Nevertheless 
direct matrix multiplication can be quite efficient if N is not too large or 
if a highly parallel computer is used. 

It is also possible to evaluate derivatives using (2.10) - (2.11). 

Indeed, (2.11) gives 


,n c 

d P N f 


dx 


n 


oo = I ( ik ) 

J |k | <N 


n. 


ikx , 




(2.16) 


where a k is given by (2.10). In this approach, a k is first evaluated by 
(2.10) and then derivatives at Xj are evaluated by (2.16). If N is a 
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highly composite integer, the two discrete Fourier transforms (2.10) and 
(2.16) can be efficiently evaluated by the FFT algorithm in 0(N log N) 
operations. Thus, evaluation of derivatives requires just two FFTs together 
with the complex multiplication by (ik) n in (2.16). 


3- POLYNOMIAL METHODS 

We now consider equation (1.1) in a finite interval -1 < x < 1 . Since 
the problem is not periodic, Fourier expansions do not yield high order 
approximations. Instead, it is preferable to use orthogonal polynomials. We 
thus take as t * ie ^asis °f % where is a polynomial of 

degree j and <}>j is zero at the appropriate boundaries. We only consider 
homogeneous boundary conditions. As before we have different spectral methods 
by choosing different {tjO and different projection operators. 

(a) Galerkin Method 

Let f(x) be a sufficiently smooth function defined in -1 < x < 1 
where f(x) vanishes at the appropriate boundaries which yields a well-posed 
problem for (1.1). Define 


N 


P N f (X) = I a k * k (x) 


k=0 


(3.1) 


where the ar.'s are chosen so that 


/ t,J (x)(P N f - f )<|> j(x)dx = 0 


j = 0, • • • ,N 


(3.2) 


for some nonnegative weight 
as 


1 3u m 

^“ (x) (sr- ' 


<o(x). We write the numerical approximation (1.2) 
L U N “ f) < f , j(x)dx =0 j =0 , • ♦ • ,N (3.3) 
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where 

N 

“» -j 0 a k (t >+k (x) 

U N <0> ‘ P N V 

As before, this gives rise to a system of N+l ordinary differential 
equations for a^t). An equivalent way is to express the p.d.e. (1.1) in 
4>-space and then truncate after N terms. 

The Legendre-Galerkin method is obtained by choosing the weight 
0)(x) = 1. The Chebyshev-Galerkin method uses the weight 
cu(x) = (1 - xV 1/2 . 

(b) Pseudospectral Chebyshev Method 

In the most common pseudospectral Chebyshev method, the 
interpolation points in the interval (-1,1) are chosen to be the extrema 

x. = cos ^ (j = 0, • • • ,N) (3.4) 

of the Nth-order Chebyshev polynomials T N (x). Here the Chebyshev polynomial 
of degree N is defined by 


T^(x) = cos(N cos 1 x)< 


(3.5) 


It follows that 


V x j> = cos 


u_ 

N * 


(3.6) 


which indicates a close relation between the pseudospectral Chebyshev and the 
pseudospectral Fourier method. In order to construct the interpolant of 
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f(x) at the point x we define the polynomials 


8j <*> = 


(1-x 2 )T'(x)(-1) 3+1 


(j = 0, • • • ,N) 


with Cq - c N - 2, Cj = 1 (1 < j < U-l). It is readily verified that 


MV ' { ik- 


The Nth-degree interpolation polynomial P N f(x) to f(x) is given by 


P N = ^ f ( x J g,(x). 

j=0 3 3 


A different way of representing P N f(x) is to use the identity 

? T k<VV x) _ i1+1 

k-0 c, " ( > ■ 

k J 

ig 

N N f(x ) N T (x.)T, (x) 

1 fOOs,00 ■ I l -J- l k i k 

j=o J J 3=0 C j k=0 c k 

, N . N f (x .)T, (x ) 

- I l V*> r 1 l — 3 3 • 

k=0 k c k j=0 c 


P N f (x) = I \ T.(x), 

k=0 K K 


2 1 ? 

a k N - E - 


C k 3=0 


It should be noted that the coefficients a k in (3.10) can be evaluated using 
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the FFT. In fact, using (3.6) in (3.10) gives 


2 i N f(x ) ^ 

\ m a— - cos -r 

O k J-° Cj 


(3.11) 


The second step in getting a pseudospectral approximation is to express 
the derivatives of P N f in terms of f(x) at the collocation points Xj. 
This can be done by differentiating either (3.8) or (3.9). With (3.8) we 
obtain 


d n P„ f(x) N n 

= l f(x.)^-g.(x) 

dx 11 j=0 j dx 11 j 


(3.12) 


so that 


dIlp N f(x k ) 
dx k 


N 

= l 

j=0 


f(x j )(D n>kj 


(3.13) 


where 


(D - ) «- ’ fir 8 k (x) L 


n jk 


dx 


(3.14) 


For example 


( Vjk 


(- 1 )* 1 * 
— X, - X. 
% 1 k 


(k * j) 


(3.15) 


'VjJ 


20 - *$> 


<Voo 


2N + 1 
6 


- <V»H 


and 

D n = ( Dl )". 


(3.16) 


12 


It should be noted from the explicit formula (3.15) that the matrix is 
not antisymmetric; also D 2 is not symmetric. These facts introduce both 
theoretical and practical difficulties in the pseudospectral Chebyshev method. 

A different way to obtain an expression for the derivative of P N f is 
to differentiate (3.9) to get 


.n 


N 


dx 


- P„ f = l 
n N , L 


k=0 


a k T k n, < x >- 


where the coefficients a^ 

dP N f 

dx 


are given by (3.10). For example. 


N N 

l a k T k (x) " l b T (x), 
c=0 K K k=0 k k 


where 


b N = °> 


b N-l 2N a N 


(3.17) 


(3.18) 


and 

\ b k = b k+2 + 2(k+1)a k+l (0 < k < N-2) . (3.19) 

In evaluating the first derivative at the collocation points x^ the FFT is 
used to evaluate 


d >N f > 

j 11 

dx 


N 

l 


k=0 



T k (x) 


where 

b k 0) ’ V 0 ‘ k < ", 


and 
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T Jn) _ An) 


(n-1) 


k b k ' " b k+2 + 2 < k+i > b k;i 




2Nb 


(n-1) 

N 


(3.20) 


The set of points Xj defined in (3.3) is not the only set used with 
pseudospectral Chebyshev approximation. For hyperbolic problems, a convenient 
alternative set of collocation points is 


it j 

y j = cos n+t 


(j = 0, • • • ,N) 


(3.21) 


Two other sets of points that are sometimes used are 


(1) _ 2tt j 

Z j C0S 2N+1 


(j = 0, • • • ,N) 


(3.22) 


and 



cos 


tt( 2 j+l) 
2N+1 


(j = 0, • • • ,N) . 


(3.23) 


(c) Pseudospectral Legendre Method 

An attractive alternative to Chebyshev polynomial expansions is 
Legendre polynomial expansions. It suffices to explain how to construct a 
pseudospectral Legendre polynomial approximation to a derivative. 

Let xq = -1, Xjj = 1, and let x^(i = 1,***,N-1) be the roots of 
q'(x), where qjj(x) is the Legendre polynomial of degree N. Given the 

values of any function f(x) at the points Xj, (j = 0,»»»,N), we construct 
the interpolating polynomial 

P N f ‘ f(X j ) 8 j <X) 


(3.24) 
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where . 

i a-*' )q£(x) 

8j(X) " " a N q N (x j ) x ~ x j 

with 

a N = N(N+1). 


Therefore, 


dx 


£ P N f 


N 9* g (x ) N 

I fOc,,) 3 = l (D £ ) jk f(x k ). 


X=X j j=0 ' k 3x* 


(3.25) 


k=0 


For example 


(M 


_ VV i 


1 jk q N (x k ) Xj - x k 


(j * k) 


(I VoO 4 “n 


(3.26) 


<Vu ’ 0 


(J<«, i * N). 


The difference between the Chebyshev and Legendre methods is evident here. 
The matrix for Legendre polynomials is nearly antisymmetric, in contrast 
to the Chebyshev matrix given in (3,15). 

By the same method, we obtain 


(D ) - - 2 ,N< *3 ) 1 

( 2> 3 k W (x . - x k ) : 


<D 2>33 


l 

3 


N 


1-x 2 

J 


1 < j,k < N-l, 
j * k 

1 < j < N-l. 


(3.27) 


.-1 


This shows that = A S A , where A is a diagonal matrix and S is 


symmetric. 
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4. ERROR EQUATION 

The differential equation we wish to solve is 

= Lu. (4«1) 


The numerical approximation given by (1.2) is 


3u 


N 


3t P N LP N U N 


(4.2) 


We define the error equation as 


1ST ‘ L “N ' <P N lP N - L)U » ’ <P B L - L >V 


(4.3) 


In the finite difference literature this is frequently called the modified 
equation. We shall now give explicit formulas for the right-hand-side of 
(4.3) for several cases of interest. 

We first consider the model hyperbolic equation 


3u _ 

3 1 3x 


-1 < X < 1, u(l,t) = 0, u(x Q ) = u Q (x). 


(4.4) 


Even though this problem has constant coefficients, nevertheless 

P„ Lu XT * L u„ because of the presence of boundaries. 

N N N 




We first _ consider Chebyshev collocation at the points x ^ = cos — , 
(see (3.4)). The method described previously satisfies (4.4) at the points 
x^, ### ,x N while at x = xq we impose the boundary condition. Since u N is 

8u n 


a polynomial of degree N and 


3x 


is a polynomial of degree N-l we have 


that 
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^ U N ^ U N 

it Jir " V* >T(t) 


(4-5) 


where Q(j is a polynomial of degree at most N and Qjg(x^) = 0, j = 1,»»»,N, 
hence Q„(x) = (l+x)T'. Comparing the coefficient of T M for both sides of 


(4.5) we see that 


1 


Alternatively, comparing both sides of (4.5) at x = x n , 


r(t) = - 


U x (1 » t) 


In conclusion the error equation for the points is 


a N ^ U N 

TT " JT + ^ (14x)t n " 3 ? ' J ~2~ (14 ^ )t n 

2N 


(4.6) 


Similarly the error equation for the points (3.21) is 


9u N 3u N, a N T N+l 9u N U x (1 ’ fc) 


3 1 3x 


N+l 3t 


(N+l) 


2 T N+1 * 


(4.7) 


and the error equation for (z^ '}, (3.22) is 


3u„ 3u 


N _ N 1 r 

TT W + 2 L t n + 


(l-bc)T' 


T = a N (t) 


2u (1) 
x 

1 + 2N 


(4.8) 


We now claim that (4.8) is also the error equation for the Galerkin 


method for (4.4). The Galerkin method satisfies (4.5) where Q and x are 


chosen so that 


17 


1 Q(x)(T (x) - l) 

A = T ( t ) / ^ d = 0 


- 1 AT 7 


j = 0, • • • ,N (4.9) 


We show that if Q(x) = 


t n (1+x)t n t o ? 

■ = —n + l T. (x) then (4.9) is 
Z k-1 K 


satisfied. 

If j = 0, Tj - 1 = 0 and the result is trivial. 
For j > 1 


A = x (t ) \ f 


1 T_(T 4 - 1) N 1 T. (T. - 1) 


0 Vi j 


dx + l J 


k v j 


dx 


- 1 ATT k< - l ATT 


. 1 T A dx N 

- x(t> -i/ —^ + 1 I Vi 

Z -1 / 2 1 k-1 J 

/I - x 


N 


■ 5 T(t) - 1 + ^ ’ °- 


It thus follows that for the constant coefficient problem (4*3) that the 
Galerkin Chebyshev and the Chebyshev collocation at the point z j^ are 
identical. 

We next consider the heat equation 


u^ = u , -1 < x < 1 , 

t xx* 9 


u(-l,t) = u(l,t) = 0, u(x,0) = u Q (x). 


(4.10) 


We now collocate at the points x^ = cos , j = 1,«**,N-1 with 
boundary conditions imposed at Xq and x^. Similar to the above procedure 
we find that the error equation for the Chebyshev collocation method is 
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9u N 

+ - . (4.11) 


3t 3x 2 


2N 


For the Galerkin method 


3u N 8 U N 

_N = N +Q(x)x(t)> 

9x 


where Q(x) Is chosen to be orthogonal to j = in the 

Chebyshev Inner product. For (4.10) we choose 
n+x - ) /■_!') Jn-v') 

<f>j = T^(x) 2 2 * ^ t ^ ien f°H° ws that the error equation is 


3u 3 u 

^ — + f- — + blT" + cT 

2 1 N j N N’ 


3t 


(4.12) 


3x 


with 


b = - 


u^.t) - (-l) N U xx (-l,t) 
2N 2 


k N 


c = 


u^q.t) + (-!) u xx (-i,t) 
2(N 2 - 1) 


5. HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS 

Before discussing stability properties of various schemes it is helpful 
to review the properties of hyperbolic partial differential equations. We 
begin with the model problem 

u = u -1 < x < 1, 0 < t < T 

t x 

(5.1) 


u(x,0) = f(x), 


u(l,t) = g(t). 
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Integration by parts leads to the inequality 

Hull 2 = g 2 (t) - u 2 (-l,t), 
or 

T T 

II U 2 ( * , T ) II + / u 2 (-l,t)dt = II f 2 ll + / g 2 (t)dt. (5.2) 

0 0 

The norms in this case are L 2 norms in space and the generalization to 
Sobolev norms is immediate. We next consider the system of equations 

= A u x -1 < x < 1, (5.3) 

where A is a P x p constant matrix. Assuming the system is strongly 

hyperbolic we can diagonalize the matrix A. Hence we replace (5.3) by 

u t = D u x -1 < x < 1, (5.4) 

where D = diag(d^ ,d^ , • • • ,d ) = (d^,d^). We order the eigenvalues d. so 

that d 1 = (d 1 ,**»,d k ) is positive and d 11 = (d fc+1 , • • • ,d ) is negative. 

The appropriate initial boundary value problem is 

u = D u -1 < x < 1 

t x 

u(x,0) = f(x) 

(5.5) 

u T (l,t) = S 1 u n (l,t) + g T (t) 

."(■1,0 = S 11 uVl.t) + g II (t), 
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where u 1 is the corresponding vector of length k and u 11 is of length 
p~ k and det(S ) ^ 0, det(S ) ^ 0. We then obtain the a priori estimate 


lu 2 (* , T ) II 2 + K / [u I (-l,t)) 2 + (u II (l,t)) 2 ]dt 

0 


< llf(x)ll 2 + K / [(g^t)) 2 + (g n (t)) 2 ]dt. 


(5.6) 


for appropriate constants Kj_, K 2 which depend on the matrices D, S 1 , S 11 . 
The above estimate holds only if the boundary conditions are dissipative, 
i.e. , 


1*11 TT 

(S 1 ) D 1 S + D 11 < 0 


II * II 
(S ) D 


s 11 + 


D > 0, 


(5.7) 


in other words, if S 1 and S** are sufficiently small. When (5,7) does not 
hold we sometimes can consider a new variable v(x,t) = e u(x,t) and then 
obtain an estimate of type (5.6) for v. Hence u(x,t) satisfies the 
inequality 


aT 

e 


< 


II u( • ,T) II 
II f (x) II 2 


2 + K / T e 2at ((u I (-l,t)) 2 + u n (l,t)) 2 dt 
0 

+ K / e 2at ( (g X (t)) 2 + (g IT (t))) 2 dt 
0 


(5.8) 


for constants K^, K 2 and a > a^. Sometimes we also need to consider a norm 
with a kernel that depends on x. 

The above estimates were all obtained for L 2 (or Sobolev) norms with a 
weight of K(x) =1. On the spectral level this is appropriate for expansions 
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in Legendre polynomials. Indeed, if one uses a Legendre pseudospectral method 
one can obtain estimates similar to (5.2), (5.6) and (5.8). However, when 
using a Chebyshev pseudospectral method it is more appropriate to consider 
Sobolev norms with a weight K(x) = (1 - x^)~ ^2 . We next show that one no 
longer gets a priori estimates of the type considered until now. Instead, one 
must rely on a weaker a priori estimate. 

We again begin with the scalar equation but with homogeneous boundary 
conditions 


u=u -1 < x < 1 

t x 

(5.9) 

u(x,0) = f(x) u(l,t) = 0. 


Following Gottlieb and Orszag [17] we consider the initial condition 


U 


lx - XqI > e 

, (5.10) 

lx - x Q | < e 


i.e., <j> £ is zero everywhere except for a £ neighborhood of Xg where <J> £ 
is a triangular function with height 1. A straightforward calculation shows 
that 

|| 4> e (x,° ) = 0(e) a) = (1 - x 2 )~ 1/2 , 


but that ll<f> £ (x,-l + e)ll 2 

l h 

of u is of order e *• . 
u(x,t) has co-norm of 


1 /, 

= 0(e z ). This shows that initially the co-norm 

However, at time 1 - e, u(x,t) = f(x+t) and so 

1 /, 

order e 4 . Hence, one cannot bound u(x,t) 


uniformly in terms of the initial conditions for the Chebyshev norm with 
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weight (l-x^) For this simple case one can overcome this difficulty by 
considering an alternative norm uj^(x) = ^ ^ > see 
Gottlieb and Orszag [17]. Heuristically , this norm helps since it is no 
longer singular at x = -1. Since the differential equation (9) only allows 
left moving waves, no difficulties arise. However, this heuristic reasoning 
also demonstrates that this cure will no longer work if there is a 
nonhomogeneous boundary condition. Waves are now created at x = 1 where the 
kernel (^(x) is still singular. In particular we consider the Initial 
boundary problem 


u = u -1 < x < 1 

t X 

(5.11) 

u(x, 0) = 0, u( 1 , t ) = <J> £ (t,t 0 ) for some t Q < T. 

T 2 

By the same argument as before / g (t)dt = 0(e). However, at a time 

2 ° 2 1/ 
t = tg + e, Hu( # ,t)ll^ and liu(*,t)lH are both of order e '2 # Hence, an 

estimate of type (5.6) cannot exist In either the to or the <o^ norm. 

Heuristically any norm which has a singular weight at either +1 or -1 will 

lead to difficulties, should the weight be zero at both ends its derivative is 

nonpositive which shows that Integration by parts will not be useful. 

The same example shows that for a system of equations even with 

homogeneous boundary conditions we cannot get an a priori estimate of type 

(5.8). Consider the system 


u 


t 


u 

X 


V 


t 


-v 

X 


-1 < X < 1 


(5.12a) 

(5.12b) 
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u(x,0) = <J> e (x,0) v(x,0) = 0 

u(l,t) = av(l,t), v(-l,t) = 3u(-l,t), a3 * 0. 


The solution for u is a left moving wave until t = 1 while v remains 

zero. At time 1 + e, u(x,t) is identically zero while 

v(x,t) = <f> e (x, ~ 1 + O- At this point v = 0(e in the norm 

2-I/9 2 - Vo 

u) = (1 - x ) z or w 2 = z which are appropriate for 

(5.12b). Hence, again one cannot bound v(x,t) uniformly in terms of initial 
condition. 

We finally show that the nonhomgeneous problem (5.1) is well-posed in any 
norm with an integrable kernel when integrations are done over x-t space. 
Define 


Q(x) = / u (x, t)dt , 
0 


(5.13) 


then 


Hence 


T 

/ (u 2 ) x dt = / ( u2 ) t dt = u2 ( x » t ^ ~ f 2 (x). 


T 

f 

0 


1 IT 

Q(x) = - / u 2 (y,T)dy + / f 2 (y)dy + / g 2 (t)dt. 


x 


0 


(5.14) 


Now integrating in space with a weight K(x) 


T 1 2 1 1 2 
/ / K(x)u (x,t)dx dt = - / K(x)dx / u (y,T)dy 


0 -1 


-1 


x 

1 


1 1 2 1 T 2 

+ / K(x)dx / f (y)dy + / K(x)dx / g (t)dt, 


-1 


-1 


0 


or 


T 1 2 1 2 T 2 

/ / K(x)u (x,t)dx dt < K 1 / f (x)dx + K. / g (t)dt. 

0-1 -1 1 0 


(5.15) 
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Furthermore setting x = -1 in (5.14) we have 


r 2 2 1 2 

/ u (-l,t)dt < / f (x)dx + / g ( t)dt . (5.16) 

0 -10 

Adding (5.15) and (5.16) we have 

T 1 T 1 T 

/ / K(x)u 2 (x,t)dx dt + / u 2 (-l,t)dt < K / f 2 (x)dx + K / g 2 (t)dt, (5.17) 

0-1 0 -1 u 0 

1 

if / K(x)dx is finite. Furthermore, one usually has the inequality 
1 2 2 2 

/ f (x) dx < K , / K(x)f (x)dx. In particular, this is true for the 

-1 -1 

Chebyshev norm. Thus we have the a priori estimate 


Hull 2 + llu(-l,t)ll 2 < C[ Ilf II 2 + 

w,x,t ’ 7 t L w,x 




(5.18) 


6. STABILITY ANALYSIS OF PSEUDOSPECTRAL SCHEMES - HYPERBOLIC EQUATIONS 

In the first section we described different spectral methods. In this 
section we shall concentrate on the collocation method and refer the reader to 
[17] for the Galerkin approach. 

For periodic domains the pseudospectral Fourier method is the most 
appropriate. We consider the model problem 

u fc = a(x) u x 0 < x < 2 tt 

( 6 . 1 ) 


u(x,0) = f(x) 


u is 2ir periodic. 
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The semi-discrete Fourier is given by 


du N 2N-1 

W ’ a(x j ) J 0 <D l>jk 


( 6 . 2 ) 


'N 


is the 2N x 2N matrix given by (2.13). 


all x we have 


llujl 


dt t a(x) ^ 2 ^ D 1 U N’ U N^' 


When a(x) is positive for 


Since is antisymmetric it follows that 


d_ 2 ? _1 U N (x>t) 
dt ji 0 a(x x ) 


0 . 


It is then straightforward to obtain an error estimate of the type 

>"<t) - V^'a.O < Cd« 2 ) 1/2 ' q Ul a q , 

where 11*11 denotes the Sobolev q norm with weight a(x). 

a, q 

When a(x) changes sign in the interval (0,2tt) the situation is less 
clear. A number of numerical studies have shown that when one uses a fixed 
number of modes then the numerical solution becomes unbounded as one 
increases t. We wish to stress that this does not mean that the numerical 
method is unstable. Stability concerns itself with a fixed time interval 
0 < t < T and lets the number of modes increase. We first consider a 
concrete example (Gottlieb, Orszag, and Turkel [15]) 

u t = a(x)u 
t v ' x 


a(x) ■ a sin x + 3 cos x + Y • 


(6.3) 



26 


One can now obtain the estimate 

II U N C t ) II x < e 1/2(,a| + ,0,)t iu N (0)H 1 . (6.4) 

We see explicitly that for 0 < t < T, we have 

H^Ct)^ < ClIu^O)^, 

however as T increases C increases exponentially. We note that for this 
example with a = 1 , 3 = Y = 0, the analytic solution is 

u(x,t) = f [ 2 tan ^(e* 1 tan^)]* 

For all t the solution is bounded when f is bounded. However, for large 
t there is a steep gradient near x = 0. Hence, when one uses a fixed number 
of modes N and lets t increase, one eventually reaches a time for which 
the mesh can no longer resolve the gradient. Hence, the main problem is not 
one of stability but rather one of resolution. 

One can show that the growth with* time is not only a difficulty with the 
collocation method but the same difficulty occurs with the Galerkin method. 
Nevertheless, Galerkin methods are used for long term meteorological 
calculations using nonlinear equations where the eigenvalues change sign. 
Hence, one may conclude that in some sense (6.1) is a harder problem than the 
nonlinear Euler equations (with smooth solutions) since (6.1) develops a sharp 
gradient as t increases. 

We also note that one can stabilize the pseudospectral algorithm for 


(6.1) by rewriting (6.1) as 
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and discretizing by 


(au) + au a u 
x x_x 

u t 2 " 2 ’ 


du XT i a u 

“ y[D(au) + aDu] - X 


dt 2 l 


T # 


(6.5) 


( 6 . 6 ) 


The approximation is now clearly axisymmetric plus a lower order term and 
hence stable for all a(x). Nevertheless computations indicate that there is 
no advantage to (6.6) over (6.2). 

Returning to (6.1) Tal-Ezer [48] has also shown that the solution is 
stable when a(x) = sin 2x. Computations show that the eigenvalues have a 
bounded (independent of N) real part for a(x) = sin kx. Presently there is 
no proof that (6.1) is stable for any a(x) using the Fourier pseudospectral 
methods. Nevertheless many computations indicate that the difficulty is one 
of growth in time and lack of resolution and not stability and convergence. 

One can alleviate the growth in time by filtering the higher modes in the 
expansion of u^. Specifically one modifies the collocation method by 
considering (Kreiss and Oliger [23] and Majda, McDonough, and Osher [31]). 


where 


N 


U N " E ®k 


P,. e 


ikx 


k=-N 


p k = \ |k|-k 0 4 




Ik | < k Q 
I k | > k Q 


(6.7) 


One can consider this as a type of artificial viscosity in Fourier space which 
does not alter the spectral accuracy of the scheme. 

For the Chebyshev method we begin with the model scalar homogeneous 


equation 


28 


u = u 

t X 


u(x,0) = f(x) 


-a < x < 1 
u(l,t) = 0. 


( 6 . 8 ) 


We consider an expansion of u in terms of Chebyshev polynomials 


N 

u N (x,t) = £ a k (t)T k (x), (6.9) 


where the a^(t) are chosen so that the equation is exact at the collocation 
points* We consider two sets of collocation points, see (3.2), (3.24) 


and 


x . 
J 


cos 


12 

N 


cos 


N+l 


j = 0 , • • • ,N 
j = 0 , • • • ,N. 


Hence, we can consider to be a polynomial of degree N. 

To prove stability we first choose an appropriate norm, 
would be 


■I || 2 7T 

II U-t II = Tf 
N 0) N 



A natural norm 


When u is a polynomial of degree N or less this is equivalent to 


2 f 1 u N (x ’ t) 

Iu n"co = / dx * 


-1 /, 2 
vl - x 


However, we have already seen that even for the partial differential equation 
one cannot find an estimate of the type 


II u(x, t ) II < Cll f (x) II . 
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Hence, we shall instead consider the norm 


iiuii 


2 


OS, 


ll(l+x)ull* 


_ N 

I (1+x )u(x ,t). 
k=0 3 3 


( 6 . 10 ) 


Again when u e P , II u II ^ = / — — — ( x>t ^ x t For simplicity we shall 

1 -1 Z 2 

vl - X 

consider the collocation points y^ and simply state the results for the x... 

Since u^ is a polynomial of degree N which satisfies u at the 
points yy u N must satisfy the differential equation 


3u 

Tt 


N 


3u 

3x 


N 


t n+i^ 

N+l 


U x (1>t) 


UjjCx.O) = f(x) 


u(l,t) = 0, (6.11) 


before proceeding we state the following lemma which is an extension of one by 
Rivlin [45]. 


Lemma 


Let u = £ a, T, (x) , then 

k=0 K k 


, , - N u(x, ) 00 

lr u dx _ 1 v k _ r 

tt ^ , N ^ r ^ 

1 / n i r\ a 


- 1 /T-J 


k-0 k l-l 


2 IN 


C 0 ‘ C N ‘ 2 

= 1 otherwise 


,( 6 . 12 ) 


in particular, if u is a polynomial of degree 4N - 1 or less 


I r 1 »»<* . i ? u( *k> 

* - 1 /TT7 c * ' a 


2N* 


(6.13) 
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We now multiply (6.11) by (l+y^)u and sum over j to get 


_ . N+ 1 ~ _ N+l 

2CSIT at (1+y j )u (y r t > ■ icrar 5 l 0 (14y j ,[u (y 3 ,1 x- (< 


since the equation is exact at the points y^, j =(),•••, N while at 
N+l, l+ yj - 0. Using (6.13) we replace the sum by an integral, and so 


N 


Jrr^r l (i+y j )u 2 (y j ,t) = i / 


1 (l+x)(u ) 


2(N+1) dt 


dx 


- 1 A 2 

✓ 1 - x 


1 1 2 A 

= - If - - u — d - x — < o. 


1 (i-x)^T x 2 


(( 


or using the lemma again 


12 12 

1 d j (H 2 )u_ dx , . 1 ; _ u dx 

- 1 - 1 <1 


< 0. 


- x)/l - x 2 


(f 


Hence, 


llu„(t)ll < Ilf II . 

N ' u)^ 


If we now integrate both sides of (6.15) with respect to t we have 


T 1 u 2 dx dt 1 (1 +x)u 2 (x,T) 1 . ,2. . 

/ / = - / S dx + I (1 ± f <*> dx 


0 -1 


(1 - X) 


/Tx 2 


- 1 A-7 


- 1 /TT7 


but 


and 


1 u 2 dx 

/ <2 / 


- 1 AT7 

i 

/ (1 + X) 
-1 


1 

-1 


2 , 
U N dx 


(1 


- x)/l + X 2 


f dx 


A - X 2 


< C / 


1 e 2 , \ 

f (x) 




.14) 
j = 

.15) 

.16) 
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so 


T luT-dxdt 1 . 2 . N 

— < C / J-W- d X . 


! I 


0 - 1 AT7 


- 1 /TT7 


Hence u N satisfies the following a priori estimate. 


Theorem 

Let u e P N be a solution of (6.8) using the collocation points 
then for some constant C 

iiu. T (t)n < cl f i 

N 7 to 1 uj^ 

and (6.17) 

T 

/ II u M ( t ) II 2 dt < Cllfll 2 . 

' N v 7 0) to 


One can also prove this theorem by using the matrix representation for 
that is appropriate for the points y^. 

We also have a similar theorem for the points Xj. 


Theorem 


Let u e be a solution of (6.8) using the collocation points x j ■ 


Then 


2N A ^ 1+5C 1^ 1 2 )u + ~ 32N “ f 2a N (t) “ a N-l (t) ^ 


j=0 


(6.18) 


< in A v ~ *j> (1 -^> f S> + — 3 2 - r 11 l 2 V°> ‘ W°>A 


j=0 


The left hand side of (6.18) is indeed a norm. The sum is strictly positive 
9 9 

for nontrivial u^ unless u^(xj,t) = 0 for j = 0,««»,N-1 and 

2 

u (x N ,t) * 0 but in this case u = C(1 - x)T' and so 2a N - a N-1 


* 0. 
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7. STABILITY ANALYSIS OF PSEUDOSPECTRAL CHEBYSHEV SCHEMES 
- PARABOLIC EQUATIONS 

The analysis of the pseudospectral Fourier method is straightforward and 
is similar to the discussion of the previous section. For the Chebyshev 
method we consider the equation 

u - u (7.1a) -1 < x < 1 

t xx 

u(l,t) = u(-l,t) * 0 (7.1b) 

u(x,0) = f(x). (7.1c) 

We now only consider the collocation points Xj (3.4). The proof that the 
Chebyshev method is stable for (7.1) for these Xj is similar to that of the 
previous section. When (7.1b) is replaced by Neumann data we derive the same 
error equation as in (4.11). We then differentiate (4.11) to obtain an 
equation for v = u x with Dirichlet boundary conditions. The same energy 
estimates can be used for v. We thus have 


Theorem 

Using the Chebyshev collocation method for (7.1a-c) with the points 
Xj = cos , j = 0,***,N, we have 

_ N 9 _ N 9 

N £ V X 1 ,<:) < N E < 7 - 2 ) 

w j=0 3 j=0 J 

When (7.1b) is replaced by u x (l,t) = u x (-l,t) we have 
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7L 

N 


N 

l 


j=o 



9u ? 
N-\Z 

3x J 


2 2 

+ N Z aj(t) 



N 


l 

j=o 


(1 


x 2)(||) 2 + N 2 «2(°) f (7.3) 


where as before is the highest coefficient of the expansion of Ujj in 

terms of Tj. 


We now consider (7.1) with the boundary condition (7.1b) replaced by the 
more general case 


au(l,t) + 0 u x (l,t) = 0 
Yu(-l,t) + 6 u (-l,t) = 0. 

X 


(7.1c) 


In this case a stability proof is not yet known for the Chebyshev collocation 
method. However, one can explicitly find the eigenfunctions. Furthermore, if 
a and 6 have the same sign while Y and 6 have opposite signs then the 
eigenvalues are real, negative and distinct [10]. 


8. TIME DISCRETIZATIONS AND ITERATIVE METHODS 

In the preceding sections we have described the construction of spectral 
approximations to the spatial operator L in (1.2). In this section we 
present some of the difficulties one faces in solving the time-dependent 
equation 

9u n 

3i- - P N LP » “N • <8 ' 1) 


There are three distinct situations that arise in practice. At times one is 
only interested in the steady state version of (8.1). For such problems one 
can use a temporarily inaccurate or even inconsistent method but one wishes to 
reach the steady state quickly. When the operator L is an elliptic operator 
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one can solve the steady state equation directly using multigrid or conjugate 
gradient methods. An efficient method for hyperbolic equations will be 
presented later. 

A second possibility is that the solution to (8.1) changes in time at a 
much slower scale then it changes in space. For these types of problems one 
may use a comparatively low order scheme to discretize the time. When the 
temporal changes are on the same scale as the spatial changes, which is 
typical of many wave propagation problems, it is not useful to use a low order 
time integration scheme. 

We now present several recent developments in the construction of time 
integration formulas for the Fourier method. In one case we consider 
hyperbolic equations where we are only interested in the steady state. The 
second case occurs when the time and space dimensions are of equal Importance. 

We consider the periodic problem 


It = f x (u) + 8 y^ u ) -'A 0 < X < TT. (8.2) 

Let D^ x denote the Fourier derivative operator in the x-direction while 
D-^y represents the y Fourier derivative. Define 


u n 

H 2x U j 


n 


u j+i,k ‘ 2u j,k + u j-i,k- 


n 


(8.3) 


We solve (8.2) by a multi-step procedure. The first two steps is a modified 
Euler approximation to (8.2). We thus have 
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and 


n+ V? _ n . At /_ ,n . _ n >> 
U j,k U j,k 2 lx j,k ly 8 j,k^ 

u . . - u“ . + At(D. f n+ ^+ d, ) 

:>k j»k v lx j,k i,y 6 J ,k > 


(8.4) 


It is readily verified that this part of the scheme is unconditionally 
unstable. We thus add a correction term which is similar to one suggested by 
Lerat [27]. Thus, the final step is 


t 1 - \ 4t> 2 H 2x )(i - TT 4!> 2 .. - -• (8-5) 


n 


4 'Ay' z.y'' 1 j,k “j,k' “j,k ~j,k* 


It is clear that once a steady state is reached that u 1 ^} = u. , = u 1 ) , . It 

j >k j »k j,k 

n+ Vo n+ Vo 

follows from (8.4) that f g z = 0 and so the steady state 

solution is independent of the time step. We note that this is not true if 

the standard Lax-Wendroff finite difference formula is used instead of 

(8.4). We next show that (8.4) - (8.5) is linearly unconditionally stable if 
2 2 

2 -Jr 2 2 tt 2 3 1 3 f 

a > P (A), a > — TT- p (B) where A = and B = 7^— . We now consider 
2 2 9 u 9u 

the scalar two- dimensional equation 


u = u + u . 
t x y 


Since the problem is periodic we can take the Fourier transform of (8.4) - 

At 

(8.5): let u be the Fourier transform of u and (5,n) the Fourier 


variables, X = ^ ^ we have 


(l + a 2 X 2 sin 2 i)(l + a 2 X 2 sin 2 y)(u n+1 - u 11 ) = [iX(?+n) - •£- (?+n) 2 ]u n (8.6) 


for -ir < 5,n < it. 


Hence, the amplification matrix is 
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G(5,n) = 


(l + “ 2 sin 2 |-)(l + « 2 X 2 sin 2 “ ^(C+n 2 ) + iX(£+n) 
(l + a 2 X 2 sin 2 f)(l + a 2 X 2 sin 2 ^ ) 


.(8.7) 


By algebra it follows that |G(£,n)l < 1 if and only if 


X (5+h) 2 < a 2 (sin 2 + sin 2 y) + sin 2 sin 2 y . 


( 8 . 8 ) 


A sufficient condition for (8.8) to hold is that 


1 /rx^ 2 * S 2 + n 2 , 2 ( . 2 ? . . 2 TU 

x (?+n) < S < a [ sin y + sin y) , 


i ♦ e . , 


2 2 2 E 

5 < 2ot sin | 


for -TT < £ < IT. 


2 TT 

It is straightforward that this inequality holds when a > . 

We next consider the other extreme when the time and space variations are 
of the same order of magnitude. We again start with the scalar equation 


u. = u 
t x 


(8.9) 


u(x,0) = f (x) u is periodic. 


We solve this using a Fourier method with 2N+1 modes. Since we wish 
spectral accuracy we wish to be able to represent the waves 

e *ikt^ - -N, • • • ,0, • • • ,N. Specifically, we assume that the solution to 
(8.9) is u(x,t) = e i7rN ( x+t ). Using leapfrog in time and letting v be the 
Fourier transform to u, the spectral leapfrog approximation to (8.9) is 
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n+1 n-1 _ . . n 

7 N " V N = 2 * iAtN V N 


<8 . io) 


where for stability 


N A t < 1 /tt . 


( 8 . 11 ) 


Solving (8.10) we find that 


v N (t> - e 


ina 


where 


sin a = iriAtN, 


iNt 

while the exact solution to (8*9) satisfies v(t) = e . Comparing v^(t) 

2 3 

and v(t) Tal-Ezer [48] found that the phase error behaves like (At) N . 
This shows that if we wish to resolve the high modes we require 


(At) 2 < C(Ax) 3 , 


( 8 . 12 ) 


which is more stringent than (8.11). Hence, the accuracy requirements demand 
a much smaller At than is allowed by stability. 

In other words if we wish to resolve N modes with a leapfrog-Fourier 
method we need not 2N+1 modes but many more modes and hence the Fourier 
method is not efficient. To resolve this difficulty it is necessary to use a 
higher order accurate time integration formula. Tal-Ezer [48] has presented a 
scheme which is efficient and has spectral accuracy in time as well as 
space. The solution to 


U t = P N LP N U 


u(0) = u 


0 


(8.13) 


is given by 
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P N LP N 

u = e u 


0 * 


(8.14) 


We approximate the solution operator by 


u XT = H (P. T LP m t)u n , 
N m v N N O’ 


(8.15) 


where 


m 


H (16R) = l i K J (R) T (0), 

“ k=0 K K 


i.e., i T^(8) is a function of i0 which we call V^XiG) then 


m 


V P N LP » '>"0 ' l J k (R) \ (P H LP » t)u O- 

k=0 


where is the k-th Chebyshev polynomial, is the Bessel function of 

order k and R is larger than the difference of the largest and smallest 

eigenvalues of P N Tal-Ezer [48] has found efficient ways of 

implementing (8.15) which do not require any complex arithmetic. He also 

shows that it is more efficient to use 2N*fl modes coupled with (8.15) in 

iuNx 

time in order to resolve e rather than using more modes and leapfrog in 

time. 


9. COMPRESSIBLE FLOW 

We now consider the application of spectral methods to the Euler 
equations, i.e., compressible inviscid fluid dynamics. The main feature of 
this system of equations is that they constitute a nonlinear conservative 
system of hyperbolic equations. Hence, the solutions will frequently include 
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shocked flows. One might suspect that global methods are not suitable for 
problems with discontinuities. However, we shall see that it is still 
possible to use spectral methods for such problems. 

Majda, McDonough and Osher [31] have shown that for a linear system of 
equations, even with constant coefficients the existence of discontinuities 
reduces the accuracy of the Fourier method if one does not specially treat the 
initial data. They also show that if one truncates the initial data using a 
Galerkin procedure one retains the spectral accuracy. For equations with 
variable coefficients one must also filter the higher modes at every time 
step. These results imply that discontinuities may reduce the global accuracy 
using an L norm. Lax [25] on the other hand has argued that although the 
formal accuracy is lost nevertheless one can recover a highly accurate 
solution. This occurs since the numerical solution contains many small scale 
oscillations. Although this destroys local accuracy one can remove the 
oscillations by an appropriate post-processor and recover the spectral 
accuracy. Thus, spectral methods have high order resolution, even in the 
presence of discontinuities. It may also be shown [14] that the spectral 
method automatically satisfies the conservation form and hence the Rankine- 
Hugoniot conditions. It is also shown in [14] that one can then extend the 
theorem of Lax and Wendroff [26] and prove that if the spectral method 
converges it will converge to a weak solution of the differential equation. 
There are indications that filter (6.7) suggested by Majda, et al. [31], 
introduces an entropy condition into the numerical solution. 

At present there are three trends in dealing with shocks. The simplest 
approach is to use a finite difference artificial viscosity to reduce the 
oscillations caused by the shock. This method reduces the formal accuracy and 
gives smeared profiles of the shock. Nevertheless, reasonable results have 
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been obtained by Taylor, et al. [51] and Hussaini and Zang [20]. A different 
approach is to use shock fitting as recommended by Moretti. Salas, Zang, and 
Hussaini [47] have used this technique for a bow shock which was mapped onto 
the boundary of the domain. Other applications are discussed in [21] and 
[22]. A third possibility [14] is to truncate the high modes at each time 
step and then to locate the shock and to filter the solution on each side of 
the shock only at the final time. It is possible to use the structure of the 
spectral method to locate the shock. Thus, we compare the given numerical 
solution in Chebyshev space with a step function in order to locate the shock 
and estimate its strength. This gives a much better shock locater than is 
possible with finite difference methods. The shock locator is based on the 
fact that spectral methods give sharp discontinuities. Hence, any artificial 
viscosity that smears discontinuities will destroy the usefulness of the shock 


locator. 
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